Coulomb impurity under magnetic field in graphene: a semiclassical approach 
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We address the problem of a Coulomb impurity in graphene in the presence of a perpendicular 
uniform magnetic field. We show that the problem can be solved below the supercritical impurity 
magnitude within the WKB approximation. Without impurity the semiclassical energies correctly 
reproduce the Landau level spectrum. For a given Landau level the WKB energy depends on the 
absolute value of angular momentum in a way which is consistent with the exact diagonalization 
result. Below the supercritical impurity magnitude, the WKB solution can be expanded as a con- 
vergent series in powers of the effective fine structure constant. Relevance of our results to validity 
of the widely used Landau level projection approximation is discussed. 



I. INTRODUCTION 

Graphene, a two-dimensional honeycomb lattice of 
carbon atoms^^ at low energies can be described by 
massless Dirac fermions.?'4. This is evident in graphene's 
Landau level structure which leads to anamolous inte- 
ger quantum Hall effect (QHE) with plateaus at a xy — 
4(n + l/2)e 2 /h (n = 0, ±1, ±2, • • • )M As sample quali- 
ties improved experiments began to reveal a large num- 
ber of additional Hall plateaus^ not expected from Lan- 
dau quantization alone; the gaps associated with these 
plateaus could only be induced by electron-electron 
interactions. 6 More recently the importance of strong 
electron-electron interactions was firmly established af- 
ter fractional quantum Hall effects (FQHE) were revealed 
by transport measurements on suspended graphene sam- 
ples^ and graphene on hexagonal Boron Nitride^ (h-BN) 
substrates J£ 

Current theories of FQHE in conventional semicon- 
ducting systems rely heavily on the notion of the projec- 
tion of interactions onto a single Landau level (i.e. in gen- 
eral cases Landau level mixing is neglected). The appro- 
priateness of the Landau level projection is not trivially 
obvious for the case of graphene. In semiconducting 2D 
electron gas this is formally achieved in the limit of a large 
magnetic field B — > oo. This is justified because Coulomb 
interaction V e _ e ~ e 2 /els scales as \f~B while the single 
particle Landau level gap hoj c scales as B. So it can be 
argued there that for large magnetic fields, interactions 
between electrons in the lowest partially-filled Landau 
level cannot induce transitions to the higher Landau lev- 
els. The projection of the interactions onto the lowest 
partially-filled Landau level is not as clearly justified in 
graphene because the single particle Landau level gaps 
in graphene also scale as ~ y/~B, which is the same as the 
interaction strength. Most FQHE theories for graphene 
have nevertheless assumed that Landau level mixing is 
an inessential complication that can be ignored.— 

For this reason, it is important to study massless Dirac 
fcrmions in the presence of Coulomb interaction and a 
quantizing magnetic field, and validate Landau level pro- 
jection approximation. If Landau level projection is a 



valid approximation, effects of Landau level mixing can 
be treated perturbatively. The simplest case would be 
a two body problem. For non-relativistic particles with 
Galilean invariance, a two body problem is equivalent 
to a one-body problem once we separate the center-of- 
mass and relative motions. However such a separation 
is not possible for Dirac fcrmions, and the two body 
problem cannot be solved analytically. In the absence 
of magnetic field solutions of two-body problems with 
zero center of mass momentum are possible^ but these 
solutions do not generalize to the present case with mag- 
netic field. Therefore as a first step we study instead 
a massless Dirac fermion in the presence of a Coulomb 
impurity and a uniform magnetic field in this work, and 
address the following question: Can Landau level mixing 
effects induced by the Coulomb impurity be treated per- 
turbatively or not? In addition to the interest in its own 
right, we note many features in this one-body problem 
such as non-linear screening and supercritical instabili- 
ties have direct generalizations in the many body case, 
for example exciton condensation or spontaneous mass 
generation^ 

Unlike the Coulomb impurity proble m 13 i 14 in the zero 
magnetic field case, this problem can not be solved ex- 
actly in closed analytic form when a uniform magnetic 
field is present. We instead apply the semi-classical WKB 
method to solve this problem. While approximate, the 
WKB method is non-perturbative in the potential. As 
we will show, it gives rise to the exact Landau level spec- 
trum in the absence of the Coulomb impurity, and nu- 
merically very accurate energy spectra in its presence 
for most cases; the latter is established by comparing 
with exact diagonalization calculation using a truncated 
Hilbert space that keeps a very large number of Landau 
levels. By expanding the WKB solutions in power series 
of Coulomb impurity strength, we show that the series 
is convergent as long as the Coulomb impurity strength 
is below the supercritical instability critical point (to be 
discussed in more detail later) , thus establishing the per- 
turbative nature of the Coulomb potential induced Lan- 
dau level mixing effects. Our results thus lend support to 
the Landau level projection approximation in this limited 
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parameter range. 

We would like to stress the importance of supercritical 
instabilities in the Coulomb impurity problem in a mag- 
netic field. Without the magnetic field, supercritical in- 
stabilities have been investigated by many authorsJ^r— 
For massless Dirac fcrmions, it is accompanied with an 
infinite number of quasi-localized resonances in the hole 
sector. When the magnetic field is added, these super- 
critical instabilities are also present at the same critcal 
value of the impurity strength g c . This is because su- 
percritical instabilities are only determined by the short 
range behavior of the effective potential. The contribu- 
tion to the effective potential induced by the presence of 
a quantizing magnetic field vanishes as r — > 0, hence it 
does not influence the short-distance part of the effec- 
tive potential. Beyond g Cl each Landau level mixes with 
the quasistationary levels and the whole Hilbert space 
can not be truncated into a single Landau level. Below 
g c , we can use the WKB approximation to solve for the 
wavcfunctions and energy spectrum of a massless Dirac 
fermion. This method also captures the characteristic 
features of the supercritical instabilities. Below g c it has 
discrete energy level solutions, which becomes continuous 
beyond the critical point g c . This signals the breakdown 
of Landau level projection. 

In Sec. II, we outline the WKB method for 2D mass- 
less Dirac particles in a uniform magnetic field, and ob- 
tain the WKB wavcfunctions and Bohr-Sommerfeld (BS) 
quantization condition for eigenenergies. The BS condi- 
tion is compared with its counterpart of Schrodinger par- 
ticle. In Sec. Ill, the WKB results are shown for cases 
with and without Coulomb impurity. We also compare 
the semiclassical energies with energies obtained from ex- 
act diagonalization. Sec. IV addresses the convergence of 
the WKB energies when expanded in powers of Coulomb 
impurity strength. Finally, we provide a detailed deriva- 
tion (using Zwaan's method) of BS condition in the Ap- 
pendix. 



II. WKB METHOD FOR GRAPHENE 

A. Outline of the Problem 

Consider the problem of a single Coulomb impurity in a 
homogeneous magnetic field perpendicular to the plane of 
graphene. Define H^^* as the Hamiltonian of the prob- 
lem with positive (negative) Coulomb impurity at the K 
(K') point of Brillouin zone. With negative Coulomb 
impurity, close to the K point, the electron quasiparticle 
states are described by the Dirac Hamiltonian 



Pauli matrices, g = Za in which Z is the impurity charge, 
a = e 2 /(nhvp) is graphene's fine structure constant, and 
k is the effective dielectric constant. (1) does not in- 
volve inter-valley scattering, because in Fourier space 
Coulomb potential behaves like l/q and is dominated 
by small q. For conventional S1O2 substrates n 2.4, 
giving a « 0.92 which is much larger than that in QED 
(awl/137). 

We use the symmetric gauge (A x , A y ) = (B/2)(y, —x). 
Resorting to the rotational symmetry of the system, the 
eigenfunctions can be written in cylindrical coordinates 
as 
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and the radial eigenequation reads 
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where Ib = y/Hc/(eB) is the magnetic length, 
e = El b I (hup) in which E is the eigenenergy of 
the Hamiltonian (1) and e is dimensionless, 1 = 0, ±1, 
±2, ... is the orbital angular momentum quantum 
number. 

Different signs of Coulomb impurity can be related by 
the operation 
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It implies that, in a certain valley, a solution to 
the Dirac equation with energy E for positive (negative) 
Coulomb impurity, has a conjugate partner a z \§>) with 
energy — E for negative (positive) Coulomb impurity. On 
the other hand, different valleys can be related by the op- 
eration 
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Hence with the same Coulomb impurity, a solution 
to the Dirac equation with energy E in valley K (K 1 ), 
has a conjugate partner o~ x |\&) with energy E in valley 
K' (K). Therefore, it is enough to solve the problem of 
negative Coulomb impurity at the K point. 
Write 
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where vp as 10 6 m/s is the Fermi velocity, the canonical 
momentum II = — i/iV + (e/c)A includes the vector po- 
tential A corresponding to the magnetic field, o~i are the 



Eq. (3) can be written as two Schrodinger-like equations, 



-u"{r) + Ui{r)u(r) = ^^(r), 



-v"(r) + U 2 {r)v(r) = w v{r), 



(7) 
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and j = 1 — 1/2 is the total angular momentum quan- 
tum number. Although we have (seemingly decoupled) 
Schrodinger-like equations (7), u(r) and v(r) are still re- 
lated to each other. The reason is that the final wave- 
function (2) is a spinor, which is the superposition of the 
states in sublattice A and B (corresponding to u(r) and 
v(r) respectively). The ratio of the two functions u(r) 
and v(r) is determined by Eq. (3). 

Morse and Feshbachi^ classified the solutions of 
second-order ordinary differential equations by types of 
singular points of the equations. With two regular (r = 0, 
yZs/e) and one irregular (r — > oo) singular points respec- 
tively, Eq. (7)'s exact solutions can not be expressed 
in closed form in terms of known special functions. At 
short distance limit Ui^(r) — > (j 2 — g 2 — l/4)/r 2 , the 
wavefunction components have the form r 7+1 / 2 , with 
i 2 . For g > g c = the parameter 7 be- 
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comes imaginary, and the wavefunction oscillates dra- 
matically towards the center. We want to point out this 
remarkable behavior of the wavefunctions at short dis- 
tance does not depend on the existence of magnetic field, 
because magnetic field related potential term has higher 
order of r dependence than other terms of the poten- 
tials in Eqs. (8), (9), and is negligible when r — > 0. The 
above phenomenon is simply the supercritical instabil- 
ity, which is already well known in graphene Coulomb 
impurity problem^ - — For the impurity problem, such 
instability signals the breakdown of the Dirac vacuum. 
Virtual electron-hole pairs are created, with negatively 
charged electrons going to infinity while the holes are 
bound to the Coulomb center (our impurity has negative 
charge). For the same problem under magnetic field, the 
virtual electrons can not go to infinity, because the effec- 
tive potential is infinite when r — > 00. When the super- 
critical instability happens, each Laudau level mixes with 
the quasistationary levels to better shield the large im- 
purity charge, and we can not truncate the whole Hilbert 
space into one Landau level. 



WKB method is one of the basic and frequently-used 
methods to solve quantum mechanics problems without 
analytic solutions. Unlike perturbation theory, WKB 
method is not connected with the smallness of poten- 
tial and thus has wider applicability range allowing one 
to study the qualitative behavior of the system. It also 
gives implicit or even explicit solutions for the energies as 
functions of parameters of the system, through which we 
can judge if the potential can be considered as perturba- 
tion from a semiclassical view. WKB method was orig- 
inally created to approximately solve one dimensional, 
or radial part of higher dimensional Schrodingcr particle 
problems. We formalize the WKB method for 2D mass- 
less Dirac particle problem below. Coulomb potential 
and uniform magnetic field are considered for our inter- 
est, but they can be replaced by any scalar and vector 
potential for general consideration. 

Writing 
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the radial Eq. (3) becomes 



F(r) 
G(r) 



(10) 



$'(r) = -£>$(r), 



(11) 



where 



D 



M _ 

r 

he 


h r 
hg 


f he 

_(M. _ 

v r 


Ib 


r 


J _ 


eB r 
2c 1 
Zc 2 


-(-- 
-(-- 


r 

E_ 


VF 






(12) 



with the total angular momentum J = Hj. Within WKB, 
we expand the solution of Eq. (11) in the form 17 ' 19 



$(r) = e^/^H^V")^), 

n=0 



(13) 



where y(r) is a scalar function and ip^>(r) are spinor 
functions. 

For matrix D in (12), the total angular momentum 
J and energy E are two conserved physical quantities, 
which are independent of h. This may lead to some confu- 
sion because, say, J equals to Hj in quantum mechanical 
treatment of the system. However, when we make % — > 
and the theory returns to classical mechanics, quantum 
number j ~ 1/H keeping the physical quantity invariant. 
Therefore, the matrix D is independent of fi. 

Inserting (13) into (11) and equating the coefficients of 
equal powers of fi, the first two equations of this set are 



iy'{r)<p(°\r) = D<pW(r), 



(14) 
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r) +iy'(r)ipW(r) = D^ {1 \r). (15) 



iy'(r) and ip^\r) are obtained as the eigenvalues and 
eigenvectors of matrix D: 



iy'(r) = h\i(r) = ±ihp(r), 

I e n i 

p(r) = 
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(16) 



where B is any constant, gi(r) is any r dependent com- 
mon factors, and they are not important for our WKB 
results. Multiplying Eq. (15) by (p^°\r) on the left helps 
us to cancel the ip^\r) depended terms. Then Eq. (15) 
becomes 



^°)(r)/''(r) = 0. 



(19) 



<pM(r) = <pi(r)=Afi(r) 



^/ sgn(S(r))(S(r)+K 
S y/sgn(S(r))(S(r) - A 4 ) 



(17) 

where subscript i = ± represents the two eigenvalues 
and their corresponding eigenvectors, S(r) = ^ 

s = sgn{-^ - f ) • sgn{°- - 
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A is any constant, fi(r) 

is the r dependent common factor which has not been 
determined yet. For a complex number z in this 
paper, we choose arg(z) in the region (— tt,tt], and 
arg(z 1 ^ 2 )=&rg(z)/2. For this reason, sgn(S(r)) inside 
the square roots of Eq. (17) can not be factored out in 
order to keep the phase difference of the wavefunctions 
in two sublattices. Since matrix D is not symmetric, left 
eigenvector fi(°\r) satisfying tp^ (r)iy' (r) = ip(°\r)D is 
introduced: 



¥>(°>(r) 



<p~i{r) = B gi {r)x 



( ^ S gn(S(r))(S(r)+X t ), -s^/ S gn(S(r))(S(r) - A,) ) , 

(18), 



Substituting (17), (18) into Eq. (19), we obtain 



/i(r) = A i (r)"5. 



(20) 



In WKB approximation, we only keep functions y(r) and 
(p(°>(r) in (13). The WKB approximate solution of Eq. 
(11) is obtained as 
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where C is a constant. The general solution $(r) could 
be written as the linear combination of $+(r) and 3>_(r): 



$(r) = ci$ + + c 2 <l>- = cip 2 e 



-i f r pdr 



V ' sgn(S)(S - ip) 
s^/sgn(S)(S + ip) 



+ c 2 p 2 e 
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y/ sgn(S)(S + ip) 
Sy/sgn(S)(S - ip) 



(22) 
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where r is redefined as a dimensionless num- 
ber representing the ratio of the distance from 
origin to magnetic lengt h Ib, also redefine 
S(r) = j It - r/2, p(r) = yftk - JJrf - (j/r - r/2) 2 
and s = sgn(e — g/r) ■ sgn(j/r — r/2) using the new 
dimensionless r, are constants fixed by boundary 
condition and normalization. 

To further obtain the BS condition for eigenenergics, 
we need to distinguish between classically allowed and 
forbidden regions. Defining WKB effective potential and 
WKB effective energy 

U efJ = (f - g 2 )/r 2 + 2eg/r + r 2 /4, (23) 

ee//-e 2 +.7, (24) 
the WKB wave number can be written as 

P = V € eff - U ef f. (25) 



U e ff captures quantitatively the behavior of the super- 
critical instability, which is originally reflected by the 
wavefunction limiting behavior r 7+1 / 2 at r — > when the 
Schrodinger-like equations (7) are considered above. As 
shown in Fig. 1, for < g < g c — the WKB effec- 
tive potential U e ff is positive infinite at both r — > and 
r — > oo, which allows us to use BS condition to obtain 
quantized energy levels and the WKB wavefunction van- 
ishes at r — > 0. For g > g c = \j\, U e ff is still positively 
infinite at r — > oo but negatively infinite at r — > 0, so the 
WKB wavefunction will oscillate as r — > 0. Just like what 
we can see from Eqs. (7), (8) and (9), the original Lan- 
dau level states mix with the quasistationary states near 
the origin, and get the chance to be closer to the impu- 
rity to screen the impurity charge when the supercritical 
instability happens. For the following WKB calculation, 
we will only consider the weak coupling region [g < g c ) 
and address the question: When are the electron states 
perturbatively connected to the states in a single Landau 
level? 



5 



Ueff 




— u 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 u r 

12 3 4 

FIG. 1: WKB effective potentials U eff = (j 2 - g 2 )/r 2 + 
2eg/r + r 2 /4 for j = 1/2. The solid line is for subcritical value 
g — 0.49; the dashed line is for supercritical value g — 0.7 
(g c = 0.5 when j = 1/2). The energy e is chosen to be 1.729, 
which is approximately the WKB energy of the 1st Landau 
level when g = 0.49 calculated in Sec. III. 

There are four solutions for the quartic equation 
p 2 (r) = 0: r = 6+ ^e 2 ~ 2g + 2~j, -e + ^ e 2 + 2g + 2j, 
e — \J e 2 — 2g + 2j and — e — \J e 2 + 2g + 2j . In the weak 
coupling region < g < g c = \j\, they are all real num- 
bers. Two of the four real solutions are positive while the 
other two are negative. The two negative solutions have 
no physical meaning, but we want to keep them for the 
calculation in Sec. III. We can label the four solutions 
a, 6, c, d and require a>b>0>c>d. Using Zwaan's 
method (in Appendix), BS condition is obtained as 

J\(r)dr = {n BS + ^^)K, (26) 

where ubs = 0, 1, 2, 9(j) is step function, 9(j) = 1 for 
j > and 9(j) = for j < 0. For the special case of 
riBS = 0, classically allowed region disappears; this corre- 
sponds to the zeroth Landau level and will be discussed 
in detail later. The BS condition of the same case for 
Schrodinger particle was obtained in Ref. [l8|. Besides 
the different forms of p(r), Dirac particle BS condition 
has additional terms 7r/2 — 9(j)n/2, instead of 7r/2. 

Compared to Zwaan's method used in Appendix, 
there is a more elegant wa y 19 ' 20 to deduce the Bohr- 
Sommerfeld condition (26) without connecting boundary 
conditions. Moreover, this method helps us to see the ori- 
gin of 9(j) in (26) straightforwardly. In Fig. 2, we draw 
the two branches of momenta ±p as functions of r. They 
join at the turning points a, b to form a single clockwise 
closed curve C in phase space. Considering any term 
of the WKB wavefunction (22), say $_(r), it must be 
single- valued after a full cycle along C. First consider the 
exponent containing J pdr. In one complete cycle, the 
phase change is § c p{r)dr. There are additional phases 
introduced at each turning point by the amplitude fac- 
tor p~ 2 . At each turning point, p will change sign, which 
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FIG. 2: Schematic diagram of the two branches of momenta 
±p in phase space. The two branches merge continuously at 
the turning points a, b. 

equals to adding a phase 7r since exp(in) = — 1. Therefore 
$_(r) gets an additional phase of — 7r/2 at each turning 
point. Another kind of additional phases is introduced at 
each point where S(r) — by the spinor factor of $_(r). 
At each such point sgn(S(r))ip(r) inside the square root 
changes sign. Each S(r) =0 point gives an additional 
phase of tt/2. Overall, the single- valuedness of the wave- 
function demands 

f 7T 7T 

j) p(r)dr - n - + k - = n BS 2n, (27) 

where /i is the number of turning points and k is the num- 
ber of points where S(r) — in one complete cycle. For 
negative j, S(r) = j jr — r/2 is always negative. For pos- 
itive j, S(r) is monotonically decreasing function, which 
equals to zero at one point in classically allowed region. 
Overall, for our case, [i — 2 and K = 29(j), so (27) returns 
to BS condition (26) directly. 

III. WKB RESULTS 

A. WKB approximation without Coulomb 
impurity and the zeroth Landau level 

In this subsection, we first turn off the Coulomb poten- 
tial, and consider the problem of one 2D Dirac particle 
in a perpendicular constant magnetic field. With exact 
solutions available, this problem enables us to compare 
WKB energies to exact solutions, and find the relation- 
ship between WKB quantum ubs a- n d energy quantum 
n. In the units we have chosen, the exact eigenenergies 
of this problem are e = ±v2n where n = 0, 1, 2, 3, . . ., 
and / > — n + 1. 

In the Bohr-Sommcrfcld quantization condition 
(26), p(r) = \/ e 2 — (j/r — r/2) 2 without Coulomb im- 
purity and the integral can be carried out as 
J b a p(r)dr - (e 2 + j>/2 - \j\n/2. For j > 0, Eq. (26) 
gives e = ±v / 2nss; thus the WKB energies are identi- 
cal to the exact energies! This also tells us nss = n 
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is the Landau level index. For j < 0, Eq. (26) gives 
e = ±^/2(ubs ~ I + 1)) which reproduces the Landau 
level spectrum when ngs = n + I — 1. Therefore, semi- 
classical energy correctly reproduces the Landau level 
spectrum. 

From the argument above, we see WKB energy for the 
zcroth Landau level is obtained when ubs = 0. However, 
in this case the two positive real roots of p 2 (r) = equal 
to each other and classically allowed region becomes one 
point, thus the integral in Eq. (26) is zero. Such situation 
never occurs for Schrodinger particles as the presence of 
the 1/2 shift in the BS condition. For Dirac particles 
such shift is zero for some cases. This is related to the 7r 
Berry phase associated with their cyclotron motions (see 
Appendix of Ref. @ for a discussion of this point). 



B. WKB approximation with Coulomb impurity 

When the Coulomb impurity is added (finite g), we can 
only carry out the energy calculation numerically. The 
integral in Eq. (26) can be expressed in terms of com- 
plete Legendre elliptic integrals of the first, F(x), second, 
E(x), and third, U(v,x) kindsi^i The Bohr-Sommerfeld 
quantization condition (26) then gives the transcendental 
equation 

~ (-1 3 + CZi - £Zb - m 2 I_!) = (n BS + — |^ )tt, 



with C = 4j + 4e 2 , = 8ge, m 2 = Aj 2 - 



(28) 
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E(x) 
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dr 
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(c 3 F( X ) 



+ 3c (b - c)U(u, x) + 3c(6 - c) 2 T 2 + (b - e) 3 T 3 ), 



where v = (a — b)/(a — c), x = \/ v ( c — d) J(b — d) and 
R(r) = y/-(r A -(r 2 +l3r + m 2 ). 

The semiclassical spectra as functions of g for Lan- 
dau levels from the —3rd to the 3rd and j up to 9/2 are 
shown in Fig. 3. For g = 0, the WKB solution coincides 
precisely with the exact solution, and the Landau level 
degeneracy of different angular momenta is lifted by fi- 
nite g. For fixed g, the smaller the higher energy. The 
reason is that j ± 1/2 stands for the Dirac particle's or- 
bital angular momentum in sublattice A(B). The particle 
with smaller angular momentum is closer to the impurity 
and feels stronger Coulomb potential. 

To assess the accuracy of the semiclassical spectra, we 
use exact diagonalization (ED) method to obtain the en- 
ergies of the same cases. For each total angular momen- 
tum labeled by j, the bases is from the —500th to the 
+500th Landau level's states when there is no impurity. 
We observe that WKB energy levels are quite close to 
the ED energy levels in Fig. 3 except for the zeroth Lan- 
dau level. First fix the value of impurity magnitude g. 
For smaller values of \j\ and energy quantum |n|, the de- 
viations from ED energy become larger. We make the 
assumption that ED results are accurate results, because 
the dimension of the bases is so large (1001D) and we 
are only considering the first few Landau levels. WKB 
method supposes the potential varies rather slowly in 
comparison to the de Broglie wavelength of the parti- 
cle. In a certain Landau level (fixed n), the particle with 
smaller \j\ is closer to the impurity, and feels steeper po- 
tential, so the WKB method becomes more inaccurate. 
On the other hand, with certain j, the particles with 
smaller \n\ have less kinetic energy (larger de Broglie 
wavelength), then the WKB method becomes more in- 
accurate too. For g = 0, WKB energy coincide with ED 
energy, with increasing g and fixed quantum numbers 
j and n, the difference becomes larger because poten- 
tial becomes steeper. Since classically allowed region is 
only one point for the zeroth Landau level, wavelength 
becomes infinite and semiclassical approximation is not 
able to give accurate results. 



IV. CONVERGENCE OF THE WKB 
SOLUTIONS 



3( X 2 (3-2u) + u(u-2)) 



In this section, we will consider the convergence of 



T3 = [ 4 (x 2 - v) (1 - v) 8 (x 2 - v) (l - v? ) F(x) ' semiclassical solved energy in the region < g < \ 

' ' Based on the transcendental equations (26), (28) and 

3 (x 2 (3 - 2u) + v(y - 2)) ^ 3x 2 - v (l + x 2 ) \ ^ analytic implicit function theorem^ we will argue that 

8 (x 2 - u ) 2 (l - v ) 2 2 (x 2 - v ) (1 - v) J ' ' semiclassical energy levels do converge in the region 

3u(x 2 (3-2v) + u(u-2)) 0<9<\j\- 

S^-^l-^ (X) ' Write f(g, e) = /» p(r)dr = (-l s + &i - PX Q - m 2 X_ 1 )/2. 

By analytic implicit function theorem, if we 



FIG. 3: (Color online) WKB energy and exact diagonalization (ED) energy levels from the negative 3rd to the positive 3rd 
Landau level and quantum number j from its minimum value in each Landau level up to 9/2. The lines are semiclassical 
energies. Negative j are in dashed lines and positive j are in regular lines. The circles (triangles) label the ED energies of 
positive (negative) j. For both WKB and ED energies, the spectra of \j\ = 1/2, 3/2, 5/2, 7/2, 9/2 are in red, brown, magenta, 
green, blue respectively. Energy spectrum with quantum number j has the range g £ [0, \j\). 



can prove for all the points satisfying g £ [0, \j\) 
and e ^= g/ ^2\j\, function /(g,e) is analytic and 
df(g,e)/dey^0, then it can be concluded that 
there exists an explicit function e(g) satisfying 
{(<?, <g))\g e [0, \m = {(g, e) e (0 < g < \j\,e ± g/^2\f\) 
1/(3, e) = i n BS + (1 - and e(g) is analytic in 

the region [0, \j\). t = g/^j2\j\ is a set of lines which 
make f(g,e) = (corresponding to the zeroth Landau 
level and being analytic obviously). The WKB results of 
nonzeroth Landau levels are definitely not in these lines. 



It can be easily seen from Eq. (28) that f(g, e) is 
analytic in the region (0 < g < e ^ g/y / 2\j\) 1 since 
a, b are two positive numbers while c, d are negative 
numbers. From the integral form of f(g,e), we obtain 
df{g, e)/dt = J 6 °[(e - V)/p]dr. Because e - V ^ in the 
classically allowed region, df(g, e)/de ^ 0. 



Therefore, semiclassical energy functions e(g) are ana- 
lytic, then convergent in the region < g < \j\. 



V. CONCLUSIONS 

In this paper we have used WKB approximation to 
study Coulomb impurity in the presence of a perpendic- 
ular uniform magnetic field in graphene. We find the 
solutions are smoothly (or perturbatively) connected to 
the states of isolated Landau level states when the im- 
purity strength is below the supercritical instability crit- 
ical point, thus lending support to the widely used Lan- 
dau level projection approximation when treating many- 
electron problems in the quantum Hall regime. The 
WKB solutions are quantitatively accurate, except for 
the 0th Landau level states. On the other hand Landau 
level mixing becomes a non-perturbative effect beyond 
the supercritical instability critical point, signaling the 
breakdown of Landau level projection approximation. 
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Appendix 

To obtain the Bohr-Sommerfeld quantization condi- 
tion, we first need to connect the WKB wavefunctions in 
classically allowed and forbidden regions. In this work we 
use the method named after Zwaan (see Ref. [23| for an ex- 
ample), which is very instructive and does not make use 
of the exact solution (like Airy function). In Sec. II. A, 
we defined WKB effective potential (23), energy (24) and 
write p in terms of them (25). The schematic diagram of 
U e ff and e e ff is plotted in Fig. 4. 

Recall Eq. (22), the WKB wavefunctions in classically 
forbidden regions I and III are 

^) = C9 -v/:-(^^!>). ,30) 

I 



or 



The above two forms of the WKB wavefunctions in the 
classically allowed region differ due to a different choice 
of the lower limits in the integral, which correspond to 
the two turning points. This is done in order to match 
the wavefunctions in the classically forbidden regions I 
and III. 

Before connecting the wavefunctions of different re- 
gions, we need to have a mathematical interlude. Write 

J ip + S - e"*/.'.* h^+S')dr +l c a , h . 

' 03) 

V e-V 

where the two equations are complex conjugate to each 
other, V(r) — g/r for our case. Constants c a ^ correspond 
to the two different lower limits a and b of each inte- 



where q(r) = ip(r). And the WKB wavefunction in clas- 




region I 5 region II a region III r 



FIG. 4: Schematic diagram of WKB effective potential U e ff 
and WKB effective energy e e // as functions of r. 



sically allowed region II can be written as 



I 

gral. It is easy to check that these constants should be 
imaginary so we write them as ic a ^ where c a ^ are real 
numbers. The relations between the square roots in the 
wavefunctions (29-32) to the ones in (33) are 



Vsgn(S) (ip + S) = e- lt VV^V\ \j f^f , (34) 

y/sgn(S)(-ip+S) = e i4 ^V^^Jl, (35) 

where t = [sgn(e — V) — sgn(S)]ir / 4. 

Now, we begin to use Zwaan's method to connect the 
wavefunctions in region II and III. Near r = a, we can 
make the linear approximation q = \Z\F \(r — a), where 
F a = d[(e — V(r)) 2 — S(r) 2 ]/ dr\ r=a . Write everything in 
complex plane 



II 



V s v / sgn{S){ip+ S) J \ 



y /sgn(S)(ip + S) 
iy/sgn(S)(-ip + S) 



(31) 



$j/(r) = B[(-ip) i e - 



fr pdr ( y/s gn(S)(-ip + S) 
i^J sgn(S)(ip + S) 



B'zi-ipyh* 



jr pdr I y /sgn(S)(ip + S) 
3^sgn(S)(-ip + S) 



(32) 
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(a) 



(6) 




V a , b = V(r = a, b). Let B[ = e~ iv ^C, B' 2 = C, Eq. 
(29) becomes 



*j/(r) = C( 



|e-V|, 



FIG. 5: Two different paths of WKB wavefunction passes 
from region III (classically forbidden region) to region II (clas- 
sically allowed region) in complex plane. 



«»[/> ~ + S')]dr + c a ~t a + i] 

sgn(e - f )««[/> + + S')]dr - c a + t a + f ] 



(41) 



r 

a = pe l( ^, j \Jr — adr = (cos^ 



. . 3 

+ isvn—(i 



(36) 

When region III and II 's wavefunctions are connected 
through the upper semicircle path as Fig. 5(a) 



q(r) -> ip(r), 



q(r)dr 



i / p(r)dr, 

J a 

r M-Ca-ip-^ qdr f V s 9n{S a )(-q + S) \ 
r(r) ~ Cq { s^sgn(S a )(q + S) J 



(37) 



V *V S 5«(5a)(ip+5) 



= cyj7^yi(i P )-2 e - 



— ip+S 
V 



-it a / ip+S 



(38) 

where 5 a ,fc = 5(r = a, 6). When we connect region III 
and IPs wavefunctions through the lower semicircle path 
as Fig. 5(b), 



q(r) ->• -ip(r), 



q(r)dr 



i / p(r)dr, 

J a 



(39) 



3>in(r) = Cq 



h„-Sa " dr ( V s 9 n (Sa)(-q + S) 



V Sy/sgn{S a ){q + S) 



V svWj(^+S) 



Cv/|e- I/|(-ip)-^e l ^ pdr 



„ — ii a ip+S 



ip+S 



where 



t a = [sgn(e - V a ) - sgn(S a )]n/4: 



(40) 
and 



Similarly, by connecting the wavefunctions of region I 
and II, we obtain 



$/j(r) =A(i 



V\ 



c°4f b (P - + S'))dr + c b -t b -l] 

sgn(e - f)cos[J b r (p + + S'))dr - c b + t b - f ] 



(42) 



where tb — [sgn(e — Vb) — sgn(Sb)]n / 4. 

Eq. (41), (42) then give us the Bohr-Sommcrfeld quan- 
tization condition 



r a i 



v 



1 V'S 1 

S')]dr-c a + c b +t a -t b = (n BS + -)7r. 

(43) 

where n B 5 = 0, 1,2, ...,p(r) = -iq(r) = [(e - g/r) 2 - (j/r 
-r/2) 2 } 1 / 2 . In region II where p 2 (r) = (e - V) 2 - S 2 = 
(e — V") 2 — (j/r —r/2) 2 > 0, e — V(r) can not be zero, 
so the signs of e — V(r) at point r = a and 6 are 
the same. It is also easy to see, for j < S(r) is 
always negative; for j > 0, S b is positive and is 
negative. Overall, we obtain t a — tb — 6(j)w/2, where 
9(j) = 1 for j > and 9(j) = for j < 0. Analyzing 
Eq. (33) by substitute r — a and b, we can always get 
- J b a [V'S/(e -V) + S']/2pdr -c a + c b =0. Finally, Eq. 
(43) is simplified as 



i; 



p(r)dr = (n B s + 



1 - 0(3) 



(44) 
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